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ABSTRACT 

We present the first detailed analysis of the mass and dynamical structure of a sample of 
six early-type lens galaxies, selected from the Sloan Lens ACS Survey, in the redshift range 
0.08 ^ z < 0.33. Both Hubble Space Telescope (HST)/ACS high-resolution imaging and 
VLT VIMOS integral field spectroscopy are available for these systems. The galaxies are 
modelled — under the assumptions of axial symmetry and two-integral stellar distribution 
function — by making use of the cauldron code, which self-consistently combines gravita- 
tional lensing and stellar dynamics, and is fully embedded within the framework of Bayesian 
statistics. The principal results of this study are: (i) all galaxies in the sample are well de- 
scribed by a simple axisymmetric power-law profile for the total density, with a logarithmic 
slope Y very close to isothermal {{y') - 1 .98+0.05 and an intrinsic spread close to 5 per cent) 
showing no evidence of evolution over the probed range of redshift; (ii) the axial ratio of the 
total density distribution is rounder than 0.65 and in all cases, except for a fast rotator, does 
not deviate significantly from the flattening of the intrinsic stellar distribution; (iii) the dark 
matter fraction within the effective radius has a lower limit of about 15 to 30 per cent; (iv) 
the sample galaxies are only mildly anisotropic, with |i5| < 0.16; (v) the physical distinction 
among slow and fast rotators, quantified by the u/cr ratio and the intrinsic angular momentum, 
is already present at z > 0.1. Altogether, early-type galaxies at z = 0.08 - 0.33 are found to be 
markedly smooth and almost isothermal systems, structurally and dynamically very similar to 
their nearby counterparts. This work confirms the effectiveness of the combined lensing and 
dynamics analysis as a powerful technique for the study of early-type galaxies beyond the 
local Universe. 

Key words: gravitational lensing — galaxies: elliptical and lenticular, cD — galaxies: kine- 
matics and dynamics — galaxies: structure. 



1 INTRODUCTION 

The currently favoured cosmological scenario, the so-called 
ACDM (cold dark matter) paradigm, has been remarkably success- 
ful at explaining the large scale structure of the Universe. In the 
non-linear regime, below several Mpc, however, the situation is less 
certain, and a full understanding of the galaxy formation and evo- 
lution processes remains a work in progress. 

Within the standard paradigm, massive early-type galaxies are 
thought to be the end-product of hierarchical merging of lower 
mass galaxies, and to be embedded in extended dark matter haloes 
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(e.g. Toomre 1977; White & Frenk 1991; Barnes 1992; Cole et al. 
2000). Numerical studies of merging galaxies (e.g. Naab et al. 
2006; Jesseit et al. 2007) have managed to reproduce a number of 
observational characteristics of massive ellipticals, and have made 
clear that stringent tests of galaxy formation models require a de- 
tailed and reliable description of the intrinsic physical properties 
of real early-type galaxies, such as their mass density distribution 
and orbital structure. Furthermore, knowledge of how these galaxy 
properties evolve through time would provide much needed infor- 
mation and even stronger constraints on the theoretical predictions. 

In the last decades, local early-type galaxies have been the ob- 
ject of substantial observational and modelling efforts. These stud- 
ies have employed a variety of tracers ranging from stellar kinemat- 
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Table 1. Basic data for the six SLACS lens galaxies analyzed in this work. The data are taken from Treu et al. (2006) and Koopmans et al. (2006). 



Galaxy name 


3 




0"SDSS 


Scff.V 


Mv 




Mb 


q*.2D 


!?PA,* 












(km s-') 


(kpc) 


(mag) 


(kpc) 


(mag) 




(deg) 


(kpc) 


(IO'OMq) 


SDSS J0037-0942 


0.1955 


0.6322 


265 ± 10 


8.48 ±0.11 


-23.11 


8.67 ±0.19 


-22.25 


0.76 


9.5 


4.77 


27.3 


SDSS J0216-0813 


0.3317 


0.5235 


332 ± 23 


16.95 ± 0.92 


-23.94 


17.28 ± 1.20 


-23.06 


0.85 


79.2 


5.49 


48.2 


SDSS J0912+0029 


0.1642 


0.3240 


313 ± 12 


14.66 ±0.23 


-23.41 


15.49 ±0.39 


-22.55 


0.67 


13.2 


4.55 


39.6 


SDSS J0959+0410 


0.1260 


0.5349 


212 ± 12 


4.50 ± 0.05 


-21.45 


4.75 ± 0.03 


-20.58 


0.68 


57.4 


2.25 


7.7 


SDSS J 1627-0053 


0.2076 


0.5241 


275 ± 12 


6.68 ± 0.06 


-22.68 


6.13 ±0.09 


-21.71 


0.85 


5.6 


4.11 


22.2 


SDSS J232 1-0939 


0.0819 


0.5324 


236 ±7 


7.93 ± 0.07 


-22.59 


8.47 ±0.11 


-21.72 


0.77 


126.5 


2.43 


11.7 



Notes: For each SLACS system we list: the redshifts z\ and Zs, of the lens galaxy and of the background source, respectively; the velocity dispersion o"sdss 
measured from the 3" diameter SDSS fibre; the effective radius Ri-g and absolute magnitude M determined by fitting de Vaucouleurs profiles to the V- and 
B-band ACS images; the isophotal axis ratio 5*,2d; the position angle )?pa,* of the major axis; the Einstein radius Reinsi and the total mass Meinsi enclosed, in 
projection, inside ^Einst- 



ics (e.g. Saglia, Berlin & Stiavelli 1992, Berlin et al. 1994, Franx 
et al. 1994, Carollo et al. 1995, Rix et al. 1997, Loewenstein & 
White 1999, Gerhard et al. 2001, Borriello et al. 2003, Thomas 
et al. 2007 and the SAURON collaboration: see e.g. de Zeeuw et al. 
2002, Emsellem et al. 2004, Cappellari et al. 2006) and kinematics 
of discrete tracers such as globular clusters (e.g. Mould et al. 1990; 
Cote et al. 2003) or planetary nebulae (e.g. Arnaboldi et al. 1996, 
Romanowsky et al. 2003; also in combination with the kinematics 
of stars: de Lorenzi et al. 2008) to hot X-ray gas (e.g. Fabbiano 
1989; Matsushita et al. 1998; Fukazawa et al. 2006; Humphrey 
et al. 2006), usually finding evidence for a dark matter halo com- 
ponent and for a total mass density profile close to isothermal (i.e. 
p,o, oc r"^) in the inner regions. 

On the other hand, because of the severe observational limi- 
tations, thorough studies of distant early-type galaxies (at redshift 
z <; 0.1) are still in their infancy. Traditional analyses based on 
dynamics alone are hindered by the lack of tracers at large radii 
and by the mass-anisotropy degeneracy, i.e. a change in the mass 
profile of the galaxy or in the anisotropy of the velocity dispersion 
tensor can both determine the same effect in the measured velocity 
dispersion map. Higher-order velocity moments, which potentially 
allow one to disentangle this degeneracy by providing additional 
constraints (Gerhard 1993; van der Marel & Franx 1993), can only 
be measured with sufficient accuracy in the inner parts of nearby 
galaxies with the current instruments. Fortunately, valuable addi- 
tional information on distant early-type galaxies can be provided 
by gravitational lensing (see e.g. Schneider, Ehlers & Falco 1992), 
when the galaxy happens to act as a gravitational lens with respect 
to a luminous background source at higher redshift. Strong gravi- 
tational lensing allows one to determine the total mass within the 
Einstein radius i?Eiiist in an accurate and almost model-independent 
way (Kochanek 1991) although, due to the mass-sheet or mass- 
slope degeneracies (Falco, Gorenstein & Shapiro 1985; Wucknitz 
2002), it does not permit the univocal recovery of the mass density 
profile of the lens. 

Gravitational lensing and stellar dynamics are particularly ef- 
fective when they are applied in combination to the analysis of 
distant early-type galaxies. The complementarity of the two ap- 
proaches is such that the mass-sheet and mass-anisotropy degen- 
eracies can to a large extent be disentangled and the mass profile 
of the lens galaxy can be robustly determined (see e.g. Koopmans 
& Treu 2002, 2003, Treu & Koopmans 2002, 2003, 2004, Barnabe 
& Koopmans 2007, Czoske et al. 2008, Czoske, Barnabe & Koop- 
mans 2008, van de Ven et al. 2008, Trott et al. 2008). Recently, 
the dedicated Sloan Lens ACS Survey (SLACS; Bolton et al. 2006; 



Treu et al. 2006; Koopmans et al. 2006; Gavazzi et al. 2007; Bolton 
et al. 2008; Gavazzi et al. 2008; Bolton et al. 2008; Treu et al. 2008) 
has discovered a large and homogeneous sample of 70 strong grav- 
itational lenses', for the most part early- type galaxies, at redshift 
between z - 0.05 and 0.51. Koopmans et al. (2006) have applied 
a joint analysis to the SLACS galaxies, finding an average total 
mass profile very close to isothermal with no sign of evolution 
to redshifts approaching unity (when including also systems from 
the Lenses Structure and Dynamics Survey, e.g. Treu & Koopmans 
2004). A limitation of the method used in those works is that lens- 
ing and dynamics are treated as independent problems, and all the 
kinematic constraints come from a single aperture- averaged value 
of the stellar velocity dispersion, which could potentially lead to 
biased results. Barnabe & Koopmans (2007, hereafter BK07) have 
expanded the technique for the combined lensing and dynamics 
analysis into a more general and self-consistent method, embedded 
within the framework of Bayesian statistics. Implemented as the 
CAULDRON algorithm, under the only assumptions of axial symme- 
try and two-integral stellar distribution function (DF) for the lens 
galaxy, this method makes full use of all the available data sets 
(i.e. the surface brightness distribution of the lensed source, and the 
surface brightness and kinematic maps of the lens galaxy) in order 
to recover the lens structure and properties in the most complete 
and reliable way allowed by the data. The ideal testing-ground for 
an in-depth analysis is represented by a sub-sample of 17 SLACS 
lenses which have been observed with the VIMOS integral-field 
unit (IFU) mounted on the VLT as part of a pilot and a large pro- 
gramme (ESO programmes 075.B-0226 and 177.B-0682, respec- 
tively; PI: Koopmans), obtaining detailed two-dimensional kine- 
matic maps (first and second velocity moments) in addition to the 
HST imaging data. The first joint study conducted with the caul- 
dron code of one of these systems, SDSS J2321 at z = 0.082, has 
been presented in Czoske et al. (2008, hereafter C08). 

In this paper, we extend the study of C08 to a total of 
six SLACS lenses for which kinematic data sets are now avail- 
able: SDSSJ0037, SDSSJ0216, SDSSJ0912, SDSSJ0959 and 
SDSSJ1627, in addition to the already mentioned SDSSJ2321. 
This sample was chosen to cover a range in redshift, mass and im- 
portance of rotation which is representative of the SLACS sample. 
Since the SLACS galaxies have been shown to be statistically indis- 
tinguishable from control samples in terms of any of their known 



A further 19 systems are possible gravitational lenses, but the multiple 
imaging is not secure. 
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Table 2. Observing log. The HR_Blue grism was used in programme 075. B- 
0226, the HR.Orange in programme 177.B-0682. 



Galaxy 


Wcxp 


Tcxp (s) 


Grism 


^lr=st [A] 


SDSS J0037 


33 


18315 


HR.Blue 


[3860,5175] 


SDSSJ0216 


14 


28 840 


HR_Orange 


[3875,5350] 


SDSSJ0912 


12 


6660 


HR.Blue 


[3860, 5295] 


SDSSJ0959 


5 


10300 


HR.Orange 


[4600, 6300] 


SDSSJ1627 


12 


24720 


HR_Orange 


[4200, 5940] 


SDSSJ2321 


15 


8 325 


HS31ue 


[5350, 5450] 



observables, such as size, luminosity, surface brightness (Bolton 
et al. 2008), location on the Fundamental Plane (Treu et al. 2006) 
and environment (Treu et al. 2008), we expect that the results of the 
combined lensing and dynamics analysis described in this work can 
be generalized to the massive early-type population, nicely com- 
plementing the work done by, e.g., the SAURON collaboration on 
lower redshift and lower mass early-type galaxies. Basic informa- 
tion on the six systems under study are listed in Table 1 . The VI- 
MOS and HST observations of these systems, together with a de- 
scription of the data reduction, will be detailed in a forthcoming 
paper (Czoske et al., in preparation). 

This paper is organized as follows; in Section 2 we give a brief 
overview of the available data sets. In Section 3 we recall the basic 
features of the cauldron algorithm and the adopted mass model. 
The results of the combined lensing and dynamics analysis of the 
SLACS subsample are presented in Sections 4 and 5, with the latter 
focusing on the recovered dynamical structure of the lenses. In Sec- 
tion 6 we summarize our findings and draw conclusions. Through- 
out this paper we adopt a concordance ACDM model described by 

= 0.3, flA = 0.7 and Hq = lOO/ikms"' IVIpc"' with h = 0.7, 
unless stated otherwise. 



2 OVERVIEW OF THE DATA SETS 
2.1 Spectroscopy 

Integral-field spectroscopy for seventeen lens systems was obtained 
using the integral-field unit of VIMOS on the VLT, UT3. All ob- 
servations, split in Observing Blocks (OB) of roughly one hour, 
including calibration, were done in service mode. 

Of the six systems analyzed in this paper, three were ob- 
served in the course of a normal ESO programme, a pilot, 075. B- 
0226 (PI: Koopmans). For this programme we used the HR-Blue 
grism with a resolution of cr = 0.8A (1.9A full width at half 
maximum, FWHIVI), covering an observed wavelength range of 
4000 to 6200 A. Each OB was split into three dithered exposures 
of 555 seconds each. For the large programme, we switched to the 
HR-Orange grism with mean resolution cr.i = 0.78 A, covering the 
range 5050 to 7460 A. Only one long exposure of 2060 seconds was 
obtained for each observing block; the number of observing blocks 
was sufficient to fill in on gaps in the data due to bad instrument 
fibres through pointing-offset between subsequent OBs. 

The data were reduced using the VIPGI package (Scodeggio 
et al. 2005; Zanichelli et al. 2005). For more details on the proce- 
dure and tests of the quality of the reduced data we refer to Czoske 
et al. (2008) and Czoske et al. (in preparation). 

The kinematic parameters v and cr were determined from the 
individual spectra using a direct pixel-fitting routine. Compared to 



Czoske et al. (2008), we have made a number of modifications in 
the algorithm. In particular, we now use almost the entire wave- 
length range that is available from the spectra; noisy parts at the 
blue and red ends of the spectra were cut off. Due to the vary- 
ing redshifts of the lenses, the rest-frame wavelength ranges and 
hence the spectral features that were used in the kinematic anal- 
ysis varied from lens to lens. The template used was a spectrum 
of the K2 giant HR 19, taken from the Indo-US survey (Valdes 
et al. 2004). The native resolution of the template spectrum is 1 A 
FWHIVI ((T,i = 0.42 A). The template is first smoothed to the instru- 
mental resolution of the VIMOS spectra, corrected to the rest frame 
of the lens. Due to the low signal-to-noise ratio of the spectra from 
individual spaxels, we assume the line-of-sight velocity distribution 
(LOSVD) to be described by a Gaussian. Tests show that including 
Gauss-Hermite terms and (van der Marel & Franx 1993) does 
not improve the fit in terms of per degree of freedom and does 
not yield robust results for and /?4 and consequently for the veloc- 
ity moments v and (see also Cappellari & Emsellem 2004). The 
larger wavelength range requires us to modify the linear correction 
function used in Czoske et al. (2008) by multiplicative and additive 
polynomial corrections (following Kelson et al. 2000). Extensive 
testing shows that choosing polynomial orders of five ensures good 
fits for the continuum shape without affecting the structure of the 
small scale absorption features. 

A number of spectral features that are not well reproduced by 
the stellar template are masked. This includes in particular the IVIg b 
line which is enhanced in the lens galaxy, possibly the result of an 
[ff/Fe] enhancement, as compared to the Galactic star HR 19 (Barth 
et al. 2002) and the Balmer lines which may be partially filled in 
by emission. 

2.2 Imaging 

We use Hubble Space Telescope (HST) imaging data from ACS and 
NICMOS to obtain information on the surface brightness distribu- 
tions of the lens galaxies and the gravitationally lensed background 
galaxies. 

ACS images taken through the F814W filter form the ba- 
sis of the lens modelling. For four of the systems, we have deep 
(full-orbit) imaging (program 10494, PI: Koopmans); for the re- 
maining two (SDSSJ2321, SDSSJ0037) we use single-exposure 
images from our snapshot program (10174, PI: Koopmans). Non- 
parametric elliptical B-spline models of the lens galaxies were sub- 
tracted off the images in order to obtain a clean representation of 
the source structure without contamination from the lens galaxy 
(Bolton et al. 2008). 

For the dynamical analysis the kinematic maps are weighted 
by the surface brightness of the lens galaxy. We use the red- 
dest band possible since this gives the most reliable representa- 
tion of the stellar light. For four of the six systems described here, 
NICIVIOS images taken through the F160W filter are available. 
For SDSS J2321 and SDSS J1627, we had instead to resort to the 
F814W ACS images. Since the lensed source is typically stronger 
in F814W than in F160W, we start from the B-spline model of the 
lens galaxy to which we add random Gaussian noise according to 
the variance map of the images. The images are convolved to the 
spatial resolution of the VIIVIOS data (the seeing limit of the ser- 
vice mode observations, 0.8 arcsec) and resampled to the grid of 
the kinematic data using swarp^ (Berlin 2008). 

^ http: //terapix. ia. fr/so£t/swarp 
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Figure 1. Best model lensed image reconstruction for the system SDSS 
J0037. From the top left-hand to bottom right-hand panel: reconstructed 
source model; HST/ACS data showing the lensed image after subtraction of 
the lens galaxy; lensed image reconstruction; residuals. In the panels, North 
is up and East is to the left. 
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Figure 2. Best dynamical model for the galaxy SDSS J0037. First row: ob- 
served surface brightness distribution, projected line-of-sight velocity and 
line-of-sight velocity dispersion. Second row: corresponding reconstructed 
quantities for the best model. Third row: residuals. In the panels, North is 
up and East is to the left. 



3 COMBINING LENSING AND DYNAMICS 

In this Section we recall the main features of the cauldron algo- 
rithm and describe the adopted family of galaxy models. The reader 
is referred to BK07 for a detailed description of the method. 
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Figure 3. Best model lensed image reconstruction for the galaxy SDSS 
J0216. Panels meaning as in Fig. 1. 




Figure 4. Best dynamical model for the galaxy SDSS J0216. Panels mean- 
ing as in Fig. 2. 



3.1 Overview of the cauldron algorithm 

The central premise of a self-consistent joint analysis is to adopt 
a total gravitational potential O (or, equivalently, the total density 
profile p, from which O is calculated via the Poisson equation), 
parametrized by a set // of non-linear parameters, and use it simul- 
taneously for both the gravitational lensing and the stellar dynam- 
ics modelling of the data. While different from a physical point of 
view, these two modelling problems can be formally expressed in 
an analogous way as a single set of coupled (regularized) linear 
equations. For any given choice of the non-linear parameters, the 
equations can be solved (in a direct, non-iterative manner) to obtain 
as the best solution for the chosen potential model: (i) the surface 
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Figure 5. Best model lensed image reconstruction for tlie galaxy SDSS 
J0912. Panels meaning as in Fig. 1. 
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Figure 6. Best dynamical model for the galaxy SDSS J0912. Panels mean- 
ing as in Fig. 2. 
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Figure 7. Best model lensed image reconstruction for the galaxy 
J0959. Panels meaning as in Fig. 1. 
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Figure 8. Best dynamical model for the galaxy SDSS J0959. Panels mean- 
ing as in Fig. 2. 



brightness distribution of the unlensed source, and (ii) the weights 
of the elementary stellar dynamics building blocks (e.g. orbits or 
two-integral components, TICs, Schwarzschild 1979; Verolme & 
de Zeeuw 2002). This linear optimization scheme is consistently 
embedded within the framework of Bayesian statistics. As a conse- 
quence, it is possible to objectively assess the probability of each 
model by means of the evidence merit function (also called the 
marginalized likelihood) and, therefore, to compare and rank differ- 
ent models (see MacKay 1992, 1999, 2003). In this way, by max- 
imizing the evidence, one recovers the set of non-linear parame- 
ters I] corresponding to the "best" potential (or density) model, i.e. 
that model which maximizes the joint posterior probability den- 
sity function (PDF), hence called maximum a posteriori (MAP) 



model. In the context of Bayesian statistics, the MAP model is 
the most plausible model in an Occam's razor sense, given the 
data and the adopted form of the regularization (the optimal level 
of the regularization is also set by the evidence). However, in a 
Bayesian context, the MAP parameters individually do not neces- 
sarily have to be probable: they have a high joint probability den- 
sity, but might only occupy a small volume in parameter space. 
This situation can arise if the MAP solution does not lie in the bulk 
of the posterior probability distribution (roughly parameter-space 
volume times likelihood density). This can be particularly severe if 
the PDF is non-symmetric in a high-dimensional space where vol- 
ume can dramatically increase with distance from the MAP solu- 
tion. In that case, the MAP solution for each parameter could eas- 
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Figure 9. Best model lensed image reconstruction for tlie galaxy SDSS 
J1627. Panels meaning as in Fig. 1. 





Figure 10. Best dynamical model for the galaxy SDSS J1627. Panels mean- 
ing as in Fig. 2. 



ily lie fully outside the PDFs of each individual parameter, when 
marginalized over all other parameters. We will see this happening 
in some cases, and, although somewhat counter-intuitive, it is fully 
consistent with Bayesian statistics and a peculiar aspect of statistics 
in high-dimensional spaces. 

As discussed in BK07, the method is extremely general and 
can in principle be applied to any potential shape. However, its 
current practical implementation, the cauldron algorithm, is more 
restricted in order to make it computationally efficient and ap- 
plies specifically to axially symmetric potentials, ^{R, z), and two- 
integral DFs / = f(E,L-), where E and L- are the two clas- 
sical integrals of motion, i.e., respectively, energy and angular 
momentum along the rotation axis. Under these assumptions, the 



dynamical model can be constructed by making use of the fast 
Monte Carlo numerical implementation by BK07 of the two- 
integral Schwarzschild method described by Cretton et al. (1999) 
and Verolme & de Zeeuw (2002), whose building blocks are not 
stellar orbits (as in the classical Schwarzschild method) but TICs.^ 
The weights map of the optimal TIC superposition that best repro- 
duces the observables, in a Bayesian sense, is yielded as an out- 
come of the joint analysis. 

The code is remarkably robust and its applicability is not dras- 
tically limited by these assumptions. Barnabe et al. (2008) have 
tested CAULDRON on a galaxy model (comprising a stellar component 
and a dark matter halo) resulting from a numerical N-body simu- 
lation of a merger process, i.e. a complex non-symmetrical system 
which departs significantly from the algorithm's assumptions. De- 
spite this, it is found that several important global properties of the 
system, including the total density slope and the dark matter frac- 
tion, are reliably recovered. 

3.2 The galaxy model 

Koopmans et al. (2006) have shown that a power-law model, de- 
spite its simplicity, seems to provide a satisfactory description for 
the total mass density profile of the inner regions of SLACS lens 
galaxies, to the level allowed by the data. This has been further 
confirmed by the COS study of SDSSJ2321, the first case where 
a fully self-consistent analysis was performed. Therefore, in the 
present work, we still adopt a power-law model. If this description 
is oversimplified for the galaxy under analysis, this will usually 
have a clearly disruptive elfect on the quality of the reconstruction, 
such as very large residuals (compared to the noise level) for the 
best model lensed image, and a patchy or strongly pixelized recon- 
structed source (Barnabe et al. 2008). 

The total mass density distribution of the galaxy is taken to be 
a power-law stratified on axisymmetric homoeoids: 



p(m) ■■ 



Po 
my' ' 



< / < 3, 



(1) 



where po is a density scale, y' will be referred to as the (logarithmic) 
slope of the density profile, and 



z 

— -I- - 

"5 



R^ jf_ 

«5 %T 



(2) 



where co and uq are length-scales and q = Co/flo- 

The (inner) gravitational potential associated with a ho- 
moeoidal density distribution p(m) is given by Chandrasekhar 
(1969) formula. In our case, for y' 2, one has 



'b{R,z) 



'-2 Jo (l+T)^Jq2 + T 

J"^" log in 
, dr, 
(1-1- T)^jq^ + T 



while for y' = 2 



(3) 



(4) 



where Oq = IrrGqUg po and 

' A TIC can be visualized as an elementary toroidal system, completely 
specified by a particular choice of energy E and axial component of the 
angular momentum L-. TICs have simple 1/i? radial density distributions 
and analytic unprojected velocity moments, and by superposing them one 
can build f(E,L-) models for arbitrary spheroidal potentials (cf. Cretton 
et al. 1999): all these characteristics contribute to make TICs particularly 
valuable and inexpensive building blocks when compared to orbits. 
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Table 3. Recovered parameters and quantities for the best power-law models of the six analyzed SLACS lens galaxies. 



Galaxy name 






y' 


g 




/dm(fie/2) 


/dm (Re) 








(deg) 








(deg) 






(10"Mo) 




SDSS J0037-0942 


65? 6 


0.434 


1.968 


0.693 


8? 8 


0.10 


0.23 


3.35 


5.40 


SDSS J0216-0813 


70?0 


0.344 


1.973 


0.816 


76? 6 


0.19 


0.17 


12.20 


9.35 


SDSS J09 12+0029 


87? 8 


0.412 


1.877 


0.672 


13?3 


0.16 


0.30 


7.40 


9.08 


SDSS J0959+0410 


80?4 


0.323 


1.873 


0.930 


64? 2 


0.25 


0.30 


0.95 


7.17 


SDSS J1627-0053 


56?4 


0.369 


2.122 


0.851 


7? 3 


0.11 


0.21 


2.23 


5.93 


SDSS J232 1-0939 


67? 8 


0.468 


2.061 


0.739 


135? 5 


0.13 


0.29 


1.98 


5.22 



Notes: We list: the four non-linear parameters, i.e. the inchnation the lens strength ao, the logarithmic slope y' and the axial ratio q; the position angle i^pa; 
the dark matter fraction /j,^ within a spherical shell of radius 0.5 and 1 (respectively), obtained under the maximum bulge hypothesis; the upper limit for 
the luminous mass Me contained inside Re ; the upper hmit for the mass-to-light ratio in the B-band. 



alii + t) al(q'- + t) 

There are three non-linear parameters in the potential to be 
determined via the evidence maximization: Oo (or equivalently, 
through equation [B4] of BK07, the adimensional lens strength q-q), 
the logarithmic slope y' and the axial ratio q. A core radius can 
be straightforwardly included in the density distribution, if neces- 
sary. Beyond these parameters, there are four additional parameters 
which determine the geometry of the observed system: the position 
angle !?pa, the inclination ; and the coordinates of the lens galaxy 
centre with respect to the sky grid. Usually, the position angle and 
the lens centre (which are very well constrained by the lensed im- 
age brightness distribution) can be accurately determined by means 
of fast preliminary explorations and kept fixed afterwards in order 
to reduce the number of free parameters during the more compu- 
tationally expensive optimization and error analysis runs. Finally, 
a proper modelling of the lensed image can occasionally require 
the introduction of two additional parameters, namely the shear 
strength and the shear angle, in order to account for external shear. 

We employ a curvature regularization (defined as in Suyu et al. 
2006 and Appendix A of BK07) for both the gravitational lensing 
and the stellar dynamics reconstructions. The level of the regular- 
ization is controlled by three so-called "hyperparameters" (one for 
lensing and two for dynamics, see discussion in BK07), whose opti- 
mal values are also set via maximization of the Bayesian evidence. 
The starting values of the hyperparameters are chosen to be quite 
large, since the convergence to the maximum is found to be faster 
when starting from an overregularized system. 



4 ANALYSIS AND RESULTS 
4.1 Best model reconstruction 

We have applied the combined cauldron code to the analysis of the 
six SLACS lens galaxies in our current sample with available kine- 
matic maps. The recovered non-linear parameters for each galaxy 
best model are presented in Table 3. The listed parameters are the 
inclination ; (expressed in degrees), the lens strength ao, the log- 
arithmic density slope y' and the axial ratio q of the total density 
distribution. The Bayesian statistical errors on the parameters (i.e. 
the posterior probability distributions) are presented in Section 4.3. 

Our analysis shows that, given the current data, there is no 
need to include external shear or core radius in the modelling of 
any of the six galaxies: there is no significant improvement in the 
evidence when these parameters are allowed to vary, and their final 



values are found to be very close to zero. As mentioned in Sec- 
tion 3, for each system the lens centre and position angle are eval- 
uated in a preliminary run and then kept constant. The best model 
position angles, relative to the total mass distribution, are found to 
depart less than 10° (and, with the exception of SDSS J0959 and 
SDSS J2321, less than 3°) from the observed values obtained from 
the light distribution. This suggests that there is at most a small 
misalignment between the dark and luminous components. A sim- 
ilar conclusion was also drawn in Koopmans et al. (2006) and used 
to set an upper limit on the level of external shear. 



The reconstructed observables corresponding to the best 
model for each galaxy (with the exception of SDSSJ2321, for 
which we refer to the plots presented in COS) are shown and com- 
pared to the data in Figs 1 to 10. The reconstruction of both lensing 
and kinematic quantities appears in general to be very accurate. The 
residuals in the reconstructed lensed image are typically consistent 
with the noise level, and there is no sign of substantial discrep- 
ancies: this indicates that the underlying total density distribution 
is not significantly more complex (e.g. strongly triaxial) than the 
adopted axisymmetric power-law model, as discussed in Barnabe 
et al. (2008). The surface brightness distribution is also well recon- 
structed. The low-level ripples which are sometimes visible in the 
residuals map are due to the discrete nature of the TICs, and can 
be easily remedied by increasing the number of the employed TIC 
components, at the cost of a considerable slow-down of the opti- 
mization process, and without changing significantly the results. In 
the case of SDSS J0037, SDSS J0216 and SDSSJ0959, the resid- 
uals clearly reveal the presence of the lensed image, which is not 
prominent in the surface brightness data map. The situation is more 
complicated for the kinematic maps, where the noise level is higher 
and the data sets composed by a smaller number of pixels (with as 
few as ten pixels for SDSS J0912). The models are fairly success- 
ful in reproducing the observables, in particular the velocity maps, 
with the exception of the velocity dispersion maps of SDSS J1627 
(the reconstructed ctioj profile declines too rapidly) and possibly 
SDSS J0959 (there appears to be a stronger gradient in the data than 
in the model). Such difficulties in reproducing the velocity disper- 
sion maps of these two systems might well reflect the shortcomings 
of the two-integral DF assumption (which implies isotropic veloc- 
ity dispersion tensor in the meridional plane, i.e. o-f, = cr-), whereas 
a more sophisticated three-integral dynamical modelling could be 
required. We discuss this issue further in Section 6. 
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Figure 11. Spherically averaged mass distributions for the sample galaxies. Thick solid line: total mass profile obtained from the best reconstructed model; 
thick dashed line: luminous mass profile (calculated under the maximum bulge hypothesis); vertical dotted hne: effective radius; vertical dashed hne: Einstein 
radius; vertical dash-dotted line: outermost boundary of the kinematic data. 



4.2 Mass distribution and dark matter fraction 

For each galaxy, we have calculated from the best model the spher- 
ically averaged total mass profile, shown in Fig. 11. In order to 
assess the dark matter content of the analyzed systems, one also 
requires the corresponding profile for the luminous component, 
which can be obtained from the best-model reconstructed stellar 
DF. However, within the framework of the method, stars are just 
tracers of the total potential, and therefore an additional assumption 
is needed in order to constrain normalization for the luminous mass 
distribution (which is arbitrary unless the stellar mass-to-light ratio 
can be determined independently). Following COS, in this work we 
adopt the so-called maximum bulge approach, which consists of 
maximizing the contribution of the luminous component. In other 
words, the stellar density profile is maximally rescaled without ex- 
ceeding the total density profile'*, positing a position-independent 
stellar mass-to-light ratio. In real galaxies the stellar mass-to-light 
ratio might not be uniform, although, based on observed colour 
gradients, the effect is not expected to be strong (e.g. Kronawitter 
et al. 2000). This procedure provides a uniform way to determine 
a lower limit for the dark matter fraction in the sample galaxies. 
Moreover, as shown in Barnabe et al. (2008), the maximum bulge 
approach is quite robust and allows a reliable determination of the 
dark matter fraction (within approximately 10 per cent of the total 



mass) even when the model assumptions are violated. The spher- 
ically averaged stellar mass profiles obtained under the maximum 
bulge hypothesis are also presented in Fig. II, where we also indi- 
cate the three-dimensional radii corresponding to the effective ra- 
dius R^, the Einstein radius i^Einst and the outermost boundary of the 
kinematic data R^^^. The boundary of the surface brightness map is 
at least comparable to R^^,^, and often larger'. Although the spatial 
coverage of the lensing and kinematic data, in most cases, does not 
extend up to R^ (the only exception being SDSS J0959, for which 
i^itin ~ \.2R^), it should be noted that the more distant regions of 
the galaxy which are situated along the line of sight — and there- 
fore observed in projection — also contribute fairly significantly in 
constraining the mass model, as extensively discussed in COS. 

We find that the total mass profile closely follows the light 
in the very inner regions, which are presumably dominated by the 
stellar component, while dark matter typically starts playing a role 
in the vicinity of the (three-dimensional) radius r = R^jl, where its 
contribution in total mass is of order 10 to 25 per cent, and becomes 
progressively more important when moving outwards (the system 
SDSS J02I6, however, constitutes an exception, with its dark mat- 
ter fraction remaining roughly constant over the probed region for 
r > 10 kpc). Within a sphere of radius r = R^, approximately 15 
to 30 per cent of the mass is dark. This result is in general agree- 



This approach is effectively the early-type galaxies equivalent of the clas- 
sical maximum disk hypothesis (van Albada & Sancisi 1986) frequently 
used in modelHng the rotation curves of late-type galaxies. 



' The surface brightness maps considered for the combined analysis and 
presented in Figs 2, 4, 6, 8 and 10 are actually cut-outs of HST images 
extending over 10 arcsec. 
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ment with the conclusions of previous dynamical studies of early- 
type galaxies in the local Universe, in particular the analysis of 
21 nearly round and slowly-rotating ellipticals by Gerhard et al. 
(2001), the modelling of 25 SAURON systems (under the assump- 
tion that mass follows light, Cappellari et al. 2006), and the study 
by Thomas et al. (2007) of 17 early-type galaxies in the Coma clus- 
ter. 

From the value Me of the luminous mass inside the effective 
radius, obtained under the maximum bulge hypothesis, we also cal- 
culate for each system the corresponding upper limit for the stellar 
mass-to-light ratio (see Table 3), finding 5 < (M/L)q_s < 9. This 
is in agreement with stellar population studies, e.g. Trujillo et al. 
(2004). 

4.3 Error analysis 

In this Section we present, for each galaxy in the sample, the corre- 
sponding model uncertainties, i.e. the errors on the recovered non- 
linear parameters i, q-q, y' and q. The uncertainties are calculated 
within the framework of Bayesian statistics by making use of the 
recently developed nested sampling technique (Skilling 2004, Sivia 
& Skilling 2006; see also Vegetti & Koopmans 2008 for the first 
astrophysical application in the context of gravitational lensing). 
Nested sampling is a Monte Carlo method aimed at calculating the 
Bayesian evidence, i.e. the fundamental quantity for model compar- 
ison, in a computationally efficient way. The marginalized posterior 
probability distribution functions (PDFs) of the model parameters, 
which are used to estimate the uncertainties, are obtained as very 
valuable by-products of the method. 

Within the context of Bayesian statistics, a priori assumptions 
or knowledge on each parameter ri, are made explicit and formal- 
ized by defining the prior function p(T]i). We assign a uniform prior 
within the interval 26rii, symmetrical around the recovered best 
model value rif,j and wide enough to include the bulk of the likeli- 
hood (very conservative estimates of if/yi are obtained by means of 
fast preliminary runs), that is: 



constant for |;7b.i - jyil < Srj, 



Pin.) ■ 



(6) 







for \rib.,-r],\>Sr],. 



This choice of an uniform prior is aimed at formalizing the ab- 
sence of any a priori information within the interval 2 Srji (see e.g. 
Cousins 1995). We find, however, that the errors on the parameters 
are very small in comparison with drji, so that the prior, largely in- 
dependently of the adopted functional form, is nearly constant over 
the likelihood. Therefore, the specific choice for p(rii) is not critical 
in our case. 

For each analyzed galaxy, the histograms in Fig. 12 show the 
marginalized posterior PDF of the power-law model non-linear pa- 
rameters. Because of the marginalization involved in their evalua- 
tion, these distributions constitute the most conservative estimate 
of statistical errors on the parameters, given the data and all the 
assumptions (i.e. positing that the adopted model is the "true" de- 
scription underlying the data; cf. MacKay 1992). These errors are 
relatively small, as a consequence of the numerous constraints pro- 
vided by the data: the typical data set for each of the sample galax- 
ies consists of ~ 10"* data points or more, most of them (in the 
lensed image and surface brightness maps) characterized by fairly 
good signal to noise level. Furthermore, the maximum z/i of the 
posterior PDF (and, more generally, the bulk of the posterior proba- 
bility) for the ;-th parameter is often found to be somewhat skewed 
with respect to the corresponding best model value r], \,. This is a 



well-known projection effect arising from the marginalization over 
a single parameter of a complicated high-dimensional multivariate 
function such as the total posterior PDF, as discussed above (§ 3.1). 

The analysis conducted in this Section does not take into ac- 
count systematic uncertainties, which are frequently larger than the 
statistical errors, and more difficult to quantify. They can arise from 
a variety of sources, including incorrect modelling assumptions and 
problems associated with the generation of the data sets (see e.g. 
Marshall et al. 2007 for an in-depth treatment of the systematic 
uncertainties connected with the lens galaxy subtraction process 
and the incomplete knowledge of the PSF). The study of Barnabe 
et al. (2008) provides a more quantitative feel for the systematic er- 
rors introduced by the adoption of an oversimplified galaxy model, 
showing that even in a quite extreme case (where the reconstruction 
of lensing observables is clearly unsatisfactory) the systematic er- 
ror on y' remains within about 10 per cent. Other parameters such 
as inclination and oblateness, however, are less robust and actually 
become ill-defined if the assumption of axial symmetry does not 
hold. We note, however, that in none of the six systems under study 
are the model residuals as severe as in the simulations in Barnabe 
et al. (2008). Therefore, we expect systematic uncertainties to re- 
main within a few per cent level. 



4.4 The density profile of the ensemble 

From the combined analysis, we have found that all the galaxies in 
the ensemble have a total density profile very close to isothermal, 
with an average logarithmic slope (y') = 1.98 ± 0.05, in agreement 
with the results of Koopmans et al. (2006). There is also no evi- 
dence of evolution of the density slope with redshift (see Fig. 13). 

No correlation is found between the logarithmic slope y' and 
the axial ratio q, the effective radius, the normalized Einstein radius 
(i.e., the ratio -Rninst/^e) and the aperture averaged velocity disper- 
sion CTsDSS- 

We now want to calculate, on the basis of the analyzed sys- 
tems, the intrinsic spread around the average slope. If we assume 
that the slope y[ of each galaxy is extracted from a parent Gaussian 
distribution of centre y'^ and variance cry, then the joint posterior 
probability for the sampling is given by 



ny'„(Ty\{y[])^p{y^,cry) Y] 



exp 



(y;-yc)- 

2(D--,+iSy'-) 



^2n(al, + Syf) 



(7) 



where p(y'^,cry) is the prior on y'^ and cry (for which we adopt a 
uniform distribution) and 6y[ are the Icr errors on y[, calculated by 
considering the region around the peak which contains 68 per cent 
of the posterior probability. From Eq. (7), the maximum likelihood 
solution for y'^ is simply the average slope (y'), while for CTy, is 
obtained from the equation 



(y[ -y'c)' 



■Sy': 



(8) 



Inserting the values for y[ and 6y[ determined in our analy- 
sis, we solve Eq. (8) to find, for our sample of six distant ellipti- 
cals, an intrinsic spread ayt = 0.092^jjQg' around the average log- 
arithmic slope, corresponding to 4.6^^ t P^r ce.n\.. A measure of the 
joint posterior of these quantities is provided by Fig. 14, where we 
plot 'P{y'^,ayi) and we draw contours corresponding to posterior 



(or likelihood in the case of flat prior) ratios PIP^y^k 
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Figure 12. Marginalized posterior probability distributions of the power-law model parameters (' (inclination), y' (logarithmic slope), no (lens strength) and 
q (axial ratio) for each of the analyzed systems, obtained from the nested sampling evidence exploration (see text). From top to bottom, the uncertainties 
correspond to galaxies SDSS J0037, SDSSJ0216, SDSS J0912, SDSS J0959, SDSS J1627 and SDSSJ2321. 
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Figure 13. The logarithmic slope of the total density profile plotted against 



redshift for the six early-type lens galaxies in the ensemble. The dashed line 
indicates the slope y' = 2, con'esponding to the isothermal profile. The etTor 
bars are calculated by considering the region of the marginalized posterior 
PDF for y' (see Fig. 12) which contains 99 per cent of the probability. 
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Figure 14. The map shows the joint posterior probability, given by Eq. (7). 
as a function of and cry . The cross marks the position of the maximum. 
The contours correspond to posterior ratios 'P/'Pm^x = e^^^'^^, with AJ-^ = 
1,4,9 (see text). 
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Figure 15. Axial ratio q of the total density distribution plotted against the 
intrinsic axial ratio of the luminous distribution. For each galaxy, qj, 
has been calculated from the con'esponding observed isophotal axis ratio 
|?*,2D. by adopting the best-model recovered inclination. The error bars are 
calculated by considering the region of the marginalized posterior PDF for 
q (see Fig. 12) which contains 99 per cent of the probability. 



4.5 Axial ratio of the density distribution 

The axial ratio q of the total density distribution is found to be al- 
ways rounder than ~ 0.65. It is interesting to compare this quantity 
with the intrinsic axial ratio q-^ of the luminous distribution, ob- 
tained by deprojecting the observed isophotal axial ratio q-^jD, i e. 

q, = ^1 - (1 - 9l2D)/sin'' ' (9) 

where we use the best model value for the inclination The results 
are illustrated in Fig. 15. For four of the galaxies in the sample, the 
flattenings coincide closely, while the total distribution is rounder 
in the case of SDSSJ1627 and SDSSJ0959. The discrepancy is 
particularly conspicuous for SDSS J0959, where q/q-,, = 1.4, while 
it is only ~ 1.1 for SDSSJ1627. Intriguingly, SDSS J0959 is the 
only clearly fast-rotating galaxy in the sample (see the velocity map 
in Fig. 8 and the strongly asymmetric DF in Fig. 16(d)), and has 
also peculiar dynamical properties when compared with the rest of 
the ensemble, as discussed in Section 5. 



with Ajp- = 1,4,9. We note that these contours are only for indi- 
cation, and formally have a proper meaning only in the case of a 
multivariate Gaussian, in which casejp^ is the usual chi-square. 

If we consider a different prior in Eq. (7), the outcome is only 
slightly modified. For instance, if we adopt p oc 1/cr.j,/, which for- 
malizes the absence of a priori information on the scale (i.e. the 
order of magnitude) of cry, we find a maximum likelihood value 
of 0.085. This points out that, despite the fact that we only have a 
handful of systems, the results are essentially driven by the data, 
with the choice of the prior playing only a minimal role. 



5 RECOVERED DYNAMICAL STRUCTURE 

As reviewed in Section 3, for any given total gravitational poten- 
tial the CAULDRON algorithm determines the best-fitting dynamical 
model by means of TICs superposition. Therefore, the best model 
of each galaxy has an associated best reconstructed map of the rel- 
ative TICs weights, which is a representation in the integral space 
(E, L-) of the corresponding (weighted) two-integral DF. The maps 
are shown in Fig. 16. In this work, consistently with COB, we em- 
ploy a library of 100 TICs for the dynamical modelling. The TIC 
grid is constructed by considering = 10 elements logarith- 
mically sampled in the circular radius R^, and, for each of them, 
A'l„ = 5 elements linearly sampled in angular momentum between 
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Figure 16. Best model reconstmction of the weighted two-integral DFs of the sample galaxies. The value of each pixel represents the relative contribution of 
the corresponding TIC to the total mass of the galaxy. 



= and L- = L-„,;„(^c)> mirrored in the negative L, values*. 
Even though this grid may appear coarse due to the Umited num- 
ber of TICs employed, it represents — for the applications described 
in this paper — an excellent compromise between computational ef- 
ficiency and quality of the observables reconstruction. Increasing 
the number of grid elements has the effect of improving the surface 
brightness reconstruction, but does not change significantly the val- 
ues of the recovered parameters. It also makes the reconstructed 
TIC weights map smoother, while preserving the main features al- 
ready visible in the corresponding coarse map. 

While the phase-space DP completely determines the structure 
and dynamics of a (collisionless) stellar system, it is not immedi- 
ately intuitive to interpret. Therefore, from the best reconstructed 
DFs we derive dynamical characterizations of the galaxies, such as 
the anisotropy parameters (§ 5.1) the (f/cr, e) diagram (§ 5.2), and 
the angular momentum (§ 5.3) which can be more easily related 
with the observations. 



* The grid in the radial coordinate corresponds to a grid in the energy 
Ec, where 



£c EE E(R^) = (i>{Rc,0) + ■ 



and the circular speed is given by: 
, '50) I 

' ORliR^.O) 



(10) 



(11) 



The choice for the circular radius also sets the corresponding maximum 
angular momentum along the z axis: L;.max(Sc) = Scfc(^c)- 



5.1 Anisotropy parameters 

The global anisotropy of an elliptical galaxy is related to the dis- 
tribution of its stellar orbits and is often considered an important 
indicator of the assembly mechanism of the system (see e.g. Burk- 
ert & Naab 2005). 

For an axisymmetric system, the global shape of the velocity 
dispersion tensor can be quantified by using the three anisotropy 
parameters (Cappellari et al. 2007; Binney & Tremaine 2008) 



and 



6=1 



where 



2n- 



2-7 ' 



(12) 
(13) 

(14) 
(15) 



denotes the total unordered kinetic energy in the coordinate direc- 
tion k, and cti, is the velocity dispersion along the direction k at any 
given location in the galaxy. For an isotropic system (e.g. the clas- 
sic isotropic rotator) the three parameters are all zero. Stellar sys- 
tems supported by a two-integral DF, as assumed in our case, are 
isotropic in the meridional plane, i.e. crj, = a: everywhere, which 
implies /} = 0. 

For each object in the sample, we have calculated the 
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anisotropy parameters by integrating up to half tlie effective ra- 
dius, which is the typical region inside which the galaxy models 
are more strongly constrained. The results are reported in Table 4. 
All the systems, with the exception of SDSS J0959, have slightly 
positive 6, i.e. are mildly anisotropic in the sense of having larger 
pressure parallel to the equatorial plane than perpendicular to it. 
This is quite similar to what Cappellari et al. (2007) and Thomas 
et al. (2008), using three-integral axisymmetric orbit-superposition 
codes, find for local ellipticals; however, their samples also dis- 
play a few galaxies with clearly higher anisotropy {6 ~ 0.4). The 
fast-rotating galaxy SDSS J0959, instead, is anisotropic in the op- 
posite sense, due to the fact that for this system cr^, < crj^ = cr^ 
over most of the density-weighted volume, which translates into a 
negative 6 parameter. This property is uncommon, although not un- 
precedented, for nearby early-type galaxies: two systems out of the 
19 analyzed by Thomas et al. (2008) have i5 < 0, while no case is 
reported from the SAURON sample. 

Since for our models /? = by construction, the two remain- 
ing anisotropy parameters are univocally related by Eq. 14 so that, 
for < S < 1, 7 is necessarily negative. The two-integral DF as- 
sumption in general does not hold for nearby ellipticals (e.g. Mer- 
rifield 1991; Gerssen et al. 1997; Thomas et al. 2008); if this is the 
case also for distant ellipticals, then the recovered y could be sig- 
nificantly in error. On the other hand, as shown in Barnabe et al. 
(2008), the global parameter 6 is more robust, and can be reliably 
recovered by cauldron (typically within ~ 15 per cent) even when 
the assumptions of two-integral DF and axial symmetry are both 
violated. 




0.2 0.4 0.6 

Figure 17. Model («/o", e) diagram for the six lens galaxies in our sam- 
ple (black circles): e^, is the intrinsic ellipticity of luminous distribution, 
and v/(T is calculated from the best model by integrating up to Re/2. We 
also show, for comparison, the corresponding quantities (corrected for in- 
clination; i)/cr measured within i?e) for the 24 nearby SAURON ellipticals 
studied in Cappellaii et al. (2007), divided in fast and slow rotators (blue 
and red points, respectively). The solid line shows the location of edge-on 
isotropic rotators, assuming o- = 0. 15 (see Binney 2005). 



5.2 The global and local b/ct 

The (v/cr, e) diagram provides a classic indicator to estimate the 
importance of rotation with respect to random motions in early- 
type galaxies (see Binney 1978). For each system, we calculate 
the "intrinsic" (i.e. inclination corrected) global quantity v/cr from 
the best reconstructed DF by integrating up to R^./! (the results are 
listed in Table 4), and we plot it against the intrinsic ellipticity of 
the luminous distribution = I - q^,. The diagram is presented in 
Fig. 17 and compared with the findings for 24 SAURON galaxies, 
corrected for inclination (Cappellari et al. 2007). There is a sharp 
dichotomy in the SLACS subsample between the obviously fast- 
rotating system SDSSJ0959 and the remaining galaxies, four of 
which have v/<r w 0.2 (two of these systems have clear character- 
istics of slow rotators, as discussed in Section 5.3). 

Whereas v/cr is a global quantity which provides information 
on the general properties of the galaxy, important insights on the 
characteristics of the system at different locations of the meridional 
plane are offered by its local analogue, i.e. the ratio (u^)/!?- between 
the mean rotation velocity around the z-axis and the mean velocity 
dispersion o"^ = (cr^ -I- crj -I- cr^)/3. We illustrate this quantity in 
Fig. 18 for the galaxies in our sample. For visualization purposes, 
the plot was produced from a weighted DF map of 900 elements 
(Ne = 30 and A'^, = 15, min'ored). In order to do this, we reopti- 
mized the best model for the dynamics hyperparameters only, de- 
termining the new optimal level of the regularization, while all the 
other parameters were kept fixed. This procedure amounts to some 
extent to an interpolation over the weighted DF map of Fig. 16 with 
the aim of obtaining a smoother distribution. 



5.3 Angular momentum 

Another robust way to characterize the global velocity structure of 
a galaxy is provided by its angular momentum content. For each 
galaxy in the ensemble we calculate the (mass-normalized) com- 
ponent of the angular momentum parallel to the axis of symmetry 
as 

J P* d^x 

where R is the radial coordinate, (v^) denotes the mean azimuthal 
stellar velocity at position x and is the spatial density of stars 
as obtained by the best reconstructed DF, i.e. = j fd?u. The re- 
sults obtained by integrating inside Rc/2 (the region most strongly 
constrained by the observations for all the systems) are reported in 
Table 4 in units of kpc km s"' . 

Whereas the dimensional parameter J- has a direct physical 
interpretation as angular momentum, it is not the most practical 
way to quantify and to compare the level of ordered rotation in 
elliptical galaxies. To this purpose, we define a more convenient 
adimensional parameter as: 

I p^R\{v^)\d\ 

A- - / I - . (17) 

J p.R^iv^f +<T^d'x 

This quantity is effectively the intrinsic equivalent of the 
observationally-defined parameter introduced by Emsellem 
et al. (2007) as an objective criterion for the kinematic classifica- 
tion of early-type galaxies. Analogously to X^, j- tends to unity for 
systems which display large-scale ordered rotation, and conversely 
it goes to zero for galaxies globally dominated by random motions. 
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Table 4. Recovered dynamical quantities for the six analyzed SLACS lens 
galaxies. 



Galaxy 


S 


y 


v/cr 


Jz 


Jz 


J0037 


0.16 


-0.37 


0.37 


112 


0.248 


J0216 


0.08 


-0.17 


0.22 


11 


0.116 


J0912 


0.07 


-0.15 


0.24 


-231 


0.229 


J0959 


-0.16 


0.27 


0.76 


-158 


0.645 


J 1627 


0.16 


-0.38 


0.23 


-69 


0.181 


J2321 


0.14 


-0.32 


0.20 


-5 


0.075 



Notes: For each galaxy we list: global anisotropy parameters 6 and y (/3 = 
by construction under the model assumption of two-integral DF); global 
u/cr ratio; angular momentum along the rotation axis J. (in units of kpc km 
s"'); dimensionless rotation parameter j-. The dynamical quantities have 
been calculated, for each system, within a cylindrical region of radius and 
height equal to Se/2. 

whereas the same galaxies might have a moderate u/cr ratio due 
to the presence of small-scale rotation patterns in the high-density 
central regions. 

We have computed j- for the SLACS subsample by integrat- 
ing inside half the effective radius, listing the results in Table 4. 
The small number of galaxies in the sample and the limited spa- 
tial coverage and quality of the kinematic data sets do not allow us 
to trace a sharp demarcation line between slow and fast rotators; 
nevertheless, there is a clear indication that SDSS J0959 belongs to 
the latter, while SDSS J0216 and SDSS J2321 are part of the first 
group, although they could not be straightforwardly identified as 
slow rotators on the basis of the (v/cr, e) diagram only. The remain- 
ing three galaxies lie somewhat in between. 



6 SUMMARY AND CONCLUSIONS 

We have conducted, for the first time, an in-depth analysis of the 
mass distribution and dynamical structure of a sample of massive 
early-type galaxies beyond the very local Universe, with a redshift 
range ofz = 0.08 -0.33. 

The examined systems are six early-type lens galaxies from 
the SLACS survey for which both HST/ ACS high-resolution imag- 
ing and VLT VIMOS integral field spectroscopy are available. 
These unique, high-quality data sets of early-type galaxies beyond 
the local Universe have enabled us to carry out a joint analy- 
sis of these systems, by combining gravitational lensing and stel- 
lar dynamics in a fully self-consistent way (using the specifically 
designed code cauldron). The method is completely embedded 
within the framework of Bayesian statistics, permitting an objec- 
tive data-driven determination of the "best model", given the ob- 
servations and our priors (formalized by the choice of the regular- 
ization). This technique makes it possible — under the assumptions 
of axial symmetry and two-integral stellar DF — to disentangle to a 
large extent several classical degeneracies and to effectively "dis- 
sect" the investigated galaxies, recovering their intrinsic structure. 
We summarize and discuss as follows the main results of this study: 

(i) The global density distribution of massive early-type galax- 
ies within approximately IR^ is remarkably well described by a 
simple axisymmetric power-law profile. Despite being very sensi- 
tive to the features of the underlying mass distribution (as shown by 
simulations of non axially symmetric systems, e.g. Barnabe et al. 
2008; see also Koopmans 2005 and Vegetti & Koopmans 2008), 



the lensed images can be reconstructed almost to the noise level by 
adopting a p oc m"'' model (not even a weak external shear is re- 
quired), indicating a surprising degree of smoothness in the mass 
structure of ellipticals. While this conclusion could already be en- 
visioned from the results of the Koopmans et al. (2006) study of 
SLACS lenses, here we have shown that such smooth models are 
also consistent with the observed surface brightness and kinematics 
maps. We suggest that this significantly smooth structure might be 
related to the formation mechanisms of early-type galaxies. 

(ii) The average logarithmic slope of the total mass density dis- 
tribution is (y') = 1.98 ± 0.05, with an intrinsic spread of 4.6^q t 
per cent. The galaxies in the sample have therefore a density profile 
consistent with isothermal, corresponding to flat rotation curves, in- 
side a range of Einstein radii of 0.3 - 0.6 R^. This is in agreement 
with the findings of previous studies of ellipticals both in the local 
Universe and up to redshift of 1 (e.g. Gerhard et al. 2001, Koop- 
mans et al. 2006, Thomas et al. 2007). 

(iii) There is no evidence for evolution of the logarithmic total 
density slope within the probed range of redshifts, although this is 
not surprising given the findings of the (non self-consistent) com- 
bined lensing and dynamics study of Koopmans et al. (2006). How- 
ever, it does show that we have systematics well under control. 

(iv) The shape of the total density distribution is fairly round, 
with an axial ratio q > 0.65, and does not differ much from the 
intrinsic axial ratio of the luminous distribution (obtained via de- 
projection of the observed isophotal ratio, by making use of the 
recovered best model value for the inclination). The only excep- 
tion is represented by the case of the lenticular galaxy SDSS J0959, 
which has a total density profile much rounder that the luminous 
one. 

(v) The lower limit for the dark matter fraction, calculated 
within spherical shells under the hypotheses of "maximum bulge" 
and position-independent stellar mass-to-light ratio, falls in the 
range 10-25 per cent (of the total mass) at half the effective ra- 
dius, and rises to 15 - 30 per cent at r = R^.. This is fully consistent 
with the results of several studies of the dark and luminous mass 
distribution in nearby ellipticals (in particular Gerhard et al. 2001; 
Cappellari et al. 2006; Thomas et al. 2007). 

(vi) The SLACS galaxies in our subsample are only mildly 
anisotropic, with the global parameter -0.16 < <5 < 0.16. Five out 
of six systems are slightly flattened (beyond the contribution of ro- 
tation) by having a larger pressure support parallel to the equatorial 
plane rather than perpendicular to it; the situation is reversed in the 
case of SDSS J0959, which shows negative S. 

(vii) From the inspection of the stellar velocity maps, the global 
and local u/cr ratio and, more decisively, the intrinsic rotation pa- 
rameter j- (directly related to the angular momentum of the stellar 
component of the galaxy), it is possible to objectively quantify the 
level of ordered rotation with respect to the random motions. Two 
of the systems, namely SDSSJ0216 and SDSSJ2321, are identi- 
fied as slow rotators, while SDSS J0959 unambiguously shows the 
characteristics of fast rotators, including the relatively low mass and 
luminosity: with Mb = -20.58 it is the only system in our sample 
to lie within the typical range of absolute magnitudes for fast ro- 
tators determined by Emsellem et al. (2007). The remaining three 
galaxies exhibit a moderate amount of rotation. 

Overall, the early-type lens galaxies analyzed in this paper, in 
the redshift range z = 0.08-0.33, are found to be very similar to the 
local massive ellipticals, in terms of structural and dynamical prop- 
erties of their inner regions. Our study shows, for the first time, that 
also the physical distinction between slow and fast rotators (origi- 
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Figure 18. Maps of the local {v^)/<t ratio between the mean rotation velocity around the z-axis and the mean velocity dispersion, plotted up to Rc/2 in the 
positive quadrant of the meiidional plane. 



nally revealed by the SAURON survey: Emsellem et al. 2007; Cap- 
pellari et al. 2007) is already in place at redshift > 0.1, although 
a larger sample is necessary in order to quantify this more pre- 
cisely to higher redshifts. Since a series of studies (Treu et al. 2006; 
Bolton et al. 2008; Treu et al. 2008) has shown that the SLACS sys- 
tems are statistically identical — in terms of observational properties 
and enviroimient — to non-lens galaxies of comparable size and lu- 
minosity, we can generalize the results of this work and conclude 
that at least the most massive elliptical galaxies did not experience 
any major evolutionary process in their global structural properties 
between redshift and 0.3. Since the look-back time is only 3.7 
Gyr, this might not be surprising. However, pushing lensing and 
dynamics techniques back in redshift and cosmic time is crucial if 
we ever wish to fully understand the structural evolution of early- 
type galaxies. In this paper a first step has been taken, using more 
sophisticated self-consistent techniques. This goes beyond what is 
possible with the use of lensing or dynamics alone. 

On the other hand, we also find differences between the ana- 
lyzed systems and nearby galaxies when we compare the respective 
global anisotropy parameters. Although the values of S recovered 
for the SLACS subsample fall within the typical range for local 
galaxies, we note that we do not find any systems with S > 0.20, 
which are instead quite common in the SAURON and Coma sam- 
ples. This might be due, however, to the very modest size of our 
sample, particularly in terms of fast-rotating objects. A much more 
drastic discrepancy is evident in the distributions of the anisotropy 
parameters y and plausibly due to the limitations of our assump- 
tion of two-integral DF, which imposes crj, = al at every location. 
Isotropy in the meridional plane is not observed, in general, for lo- 
cal ellipticals, where usually a\ > en (e.g. Gerssen et al. 1997; 



Thomas et al. 2008). If this applies also to more distant systems, 
then the anisotropy will not be correctly estimated by our method, 
and more sophisticated axisymmetric three-integral models might 
provide a better dynamical description of the galaxy. We foresee 
future developments of the cauldron algorithm in this direction, 
concurrently with the availability of improved kinematic data sets. 
As for the present analysis, however, the /(£, L.) assumption gen- 
erally appears to work satisfactorily, providing generally correct re- 
constructions of the observables. Furthermore, as tested in Barnabe 
et al. (2008), the current method is robust enough to recover in 
a reliable way several important global properties of the analyzed 
galaxy (such as the density slope, dark matter fraction, angular mo- 
mentum and i5 parameter, although not the flattening and the y and 
yS parameters) even if it is applied to a complex, non-symmetric 
system which departs significantly from the idealized assumptions 
of axisymmetry and two- or three-integral DF, and which is likely 
more extreme than the typical ellipticals under study. 



In future papers in this series, we plan to extend the current 
study by applying the combined analysis to all the 30 systems for 
which two-dimensional kinematic maps are or will become avail- 
able. This includes the entire sample of 17 SLACS early-type lens 
galaxies for which integral field spectroscopy has been obtained. 
Two-dimensional kinematic maps can also be obtained for a fur- 
ther 13 lens galaxies for which long-slit spectroscopic observations 
have been conducted with the instrument LRIS mounted on Keck-1, 
with the slit positions aligned with the major axis of the system and 
offset along the minor axis in order to mimic integral field spec- 
troscopy. 
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